The Finite Temperature Mott Transition in the Hubbard Model in Infinite Dimensions 
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We study the second order finite temperature Mott transition point in the fully frustrated Hubbard 
model at half filling, within Dynamical Mean Field Theory. Using quantum Monte Carlo simulations 
we show the existence of a finite temperature second order critical point by explicitly demonstrating 
the existence of a divergent susceptibility as well as by finding coexistence in the low temperature 
phase. We determine the location of the finite temperature Mott critical point in the (U, T) plane. 
Our study verifies and quantifies a scenario for the Mott transition proposed in earlier studies Jl] of 
this problem. 
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When the ratio of the strength of the electron elec- 
tron interaction to the bandwidth is increased, a metal 
insulator transition (MIT) occurs This phenomenon 
has been continuously and intensively studied in narrow 
bandwidth systems for several decades ||, with V2O3 be- 
ing the archetypal system. 

Providing a detailed theoretical description of this 
transition in systems which are not magnetically ordered 
is one of the most challenging problems in condensed 
matter physics. In recent years, great progress was 
achieved using the dynamical mean field theory (DMFT), 
a method which becomes exact in the limit of large lattice 
coordination. Within this framework, one can describe 
both metallic and insulating phases. See Q for a review. 
Though some aspects of this theory agree with the sce- 
narios put forth earlier by Brinkman and Rice and by 
Hubbard Q|| the DMFT approach suggested important 
qualitative modifications of these pictures. 

At zero temperature, in the DMFT of the fully frus- 
trated Hubbard model, upon increasing the interaction 
the paramagnetic metallic state is destroyed at a critical 
value denoted by U C 2- On the other hand, the insulat- 
ing solution that exists for large U disappears when U 
is decreased, at a different critical value U c i, at which 
the gap closes. Since U c i < U C 2 there is a region where 
two phases i.e., a paramagnetic metal and a paramag- 
netic insulator coexist. This region of coexistence natu- 
rally extends to finite temperatures, where it is bounded 
by two lines U C \{T) and U C 2(T). Consequently, at fi- 
nite temperatures, a first order metal to insulator transi- 
tion takes place at a value of the interaction U C (T), with 
U cl (T) < U C (T) < U c2 (T) §0. 

At temperatures above this transition line, two dif- 
ferent crossover regions were identified ||. The first 
crossover region, is the natural continuation of the U C 2 (T) 
line, and was characterized as the place where the low en- 
ergy resonance in the spectral function rapidly loses its 
intensity and the resistivity increases rapidly. The sec- 



ond crossover region can be considered to be the natural 
extension of the U c i(T), and was characterized as the 
line where the gap is comparable to the temperature and 
where a crossover to activated behavior in the resistiv- 
ity is seen. These crossover lines were further related 
H to the experimental observations of Mc Whan et al. 
and Kuwamoto et al. |9|,[l0| . Note that generically within 
DMFT, i.e., both at high and low temperatures , the 
destruction of the metallicity and the gap closure occur 
at different locations on the phase diagram. The only 
exception to this is the finite temperature second order 
critical point at which the two boundary lines U c i (T) and 
U C 2(T) meet and terminate. In this paper, we provide a 
description of this second order point at which the MIT 
occurs. 

Though substantial evidence indicated Q that the 
qualitative aspects of the finite temperature MIT sce- 
nario that is described above were genuine properties of 
the exact solution of the model. Recent work quantum 
Monte Carlo (QMC) work challenged this scenario |pl| . 
In Ref. (l^] it was claimed that, within the temperature 
range numerically accessible to the QMC technique, it 
was not possible to find any evidence for a finite temper- 
ature MIT point. It is very important to clarify this is- 
sue, since our earlier work has successfully predicted new 
experiments in V2O3 [B| and provided a natural interpre- 
tation of earlier data |||nj] . We describe below, how the 
MIT point can be found and we determine its location 
with reasonable accuracy, illustrating the power, but at 
the same time the subtleties of the QMC technique when 
applied in conjunction with the self-consistent DMFT 
equations. 

Besides quelling the doubts raised recently in Ref. pH 
there are other motivations for this work. As stressed by 
Nozieres (jl2), the Mott transition lacks an obvious or- 
der parameter, so it is not clear a priori which quantities 
should exhibit singular behavior. Secondly, the second 
order Mott end point is observed experimentally in V2O3 
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and NiSe x S2- x , an d its properties are the subject of 
ongoing experiments. Finally, it is important from the 
point of view of numerical studies of correlated electron 
systems. Since a full analytic solution to this problem 
is unavailable, the understanding of models in the limit 
of large lattice coordination requires numerical work and 
analytic approximations. This is a common situation in 
the field of correlated electron systems and the Hubbard 
model in infinite dimensions provides an excellent play- 
ground to test the merit of various analytic approxima- 
tions and numerical methods. 

We consider the Hubbard model on the Bethe lattice 
in the paramagnetic phase [fl3f with coordination d and 
hopping -J= in the d — » oo limit |Q at half filling. 

<i,j>iruj i 

The half-bandwidth is given by D = 2t and we set D = 1 
as a unit of energy. The chemical potential is set equal 
to y at half-filling. 

All the local correlation functions can be obtained from 
a Single Impurity Anderson model (SIAM) with a hy- 
bridization function 

A(iu n ) = J2^— (2) 

k 

provided that A(iu n ) obeys the self-consistency condi- 
tion 

t 2 G imp (iuj n )[A,a} = A(iu) n ). (3) 

Here a denotes the control parameters of the problem U 
and T. Gi mp (iui n )[A] is the finite temperature Matsub- 
ara /-electron Green's function of the SIAM 

Hsiam = e kci a c klJ + V k (c+ a f a + f+c k<7 ) + Un n n fl 

ka 

(4) 

To solve the problem, the single impurity Green's func- 
tion is calculated using QMC or exact diagonalization 
methods and (||) is used to check for self consistency 
This iterative process is continued until self con- 
sistency is achieved. At that point, Gi mp (iuj n ) coincides 
with Gi oc {iuj n ) 1 the local Green's function of the original 
lattice problem. At each iteration step we use quantum 
Monte Carlo simulations (l^] to calculate the impurity 
Green's function in imaginary time |]l7f . This algorithm, 
is widely considered to provide numerically exact solu- 
tions to the model (in the Monte Carlo sense). 

The parameters of the simulation are as follows: we 
typically perform 60,000 sweeps (1 sweep = 1 attempt to 
update each of the L pseudo-spins = L attempts where 
L is the number of time slices used), and when required 



upto 300,000 sweeps in order to minimize the enhanced 
fluctuations close to the critical point. We always set the 
length of the time-slice At = (3/L < 0.5, where j3 is the 
inverse temperature. 

As a criterion for the convergence of the numerical so- 
lution of the DMFT equations, we monitor the evolution 
with the iteration number of G(iu)\), the value of the 
Green's function at the first Matsubara frequency. This 
quantity is appropriate, as it is essentially the integral 
of G(t) (with r being the imaginary time), the quan- 
tity that is directly computed statistically in the QMC 
calculation. Also, G(iu)\) is the value of the frequency de- 
pendent Green's function that fluctuates most, thus, it 
sets an upper bound for the statistical error of the whole 
G(iui n ). We stop the iterations when the fluctuations in 
G(iu>i) become of the order of the statistical error of the 
QMC (which is controlled by the number of sweeps at a 
given U and T) and remain stable for at least about 30 
more iterations. 

In generic regions of the (£7, T) plane less than 10 iter- 
ations are sufficient to obtain a converged solution. How- 
ever, since we are looking for a critical point which is a bi- 
furcation point of (0) , the number of iterations needed to 
obtain convergence diverges, as we approach the critical 
point. This, is the usual phenomenon of critical slowing 
down associated with bifurcating solutions in recursive 
procedures and it forces us to increase the number of iter- 
ations, as we approach the critical point. For example, up 
to 300 iterations were necessary to check the convergence 
at parameter values close to the critical point. Special 
care was taken to assure that the solutions were indeed 
converged. To this end, we used two different initial seeds 
for the iterative procedure: G(t)[U — 0] that corresponds 
to a metallic state, and a G(r)[t = 0] that corresponds 
to an insulating one. These two initial Green's functions 
have qualitative different behavior a low frequencies, one 
with a finite density of states at the Fermi energy and 
the other with an insulating gap. If there is a single so- 
lution to the mean field equations, the algorithm should 
in both cases evolve towards the unique attractor. If two 
different solutions are allowed, one metallic like and one 
insulator like, it may be expected that the metallic seed 
would evolve towards the former and the insulating seed 
towards the latter. Using this method, when we find a 
converged unique solution to (|^) , we are certain that the 
algorithm has fully converged. 

To prove the existence of a finite temperature critical 
point we need to isolate a physical quantity that exhibits 
singular behavior. We view this MIT as a liquid gas 
transition, where the double occupancy (d) plays the role 
of the density, while U plays the role of the pressure and 
T the natural role of temperature. This analogy was put 
forward in an insightful paper by Castellani et al. fig] . As 
shown here their scenario is realized within the DMFT 
solution of the Hubbard model. With this analogy to the 
liquid gas transition in mind, we focus on the behavior 
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of the double occupation as a function of temperature 
and interaction strength. We define the susceptibility 
X = maxu[d(d)/dU] that, as we shall later show, is a 
quantity that diverges at the finite T second order critical 
point. 

In Fig. [j], we show our results for the double occupa- 
tion (d) as a function of U. Each set of data is obtained 
at a constant T, and we have checked the uniqueness of 
the solutions using the procedure described above. 

The most remarkable feature is the rapid variation of 
(d) with the interaction U at the lower temperatures [j7j . 
This variation is a direct consequence of a divergence in 
the susceptibility x defined earlier. In order to obtain 
the susceptibility from the discrete set of data at each 
temperature, we fit the numerical results with analytic 
expressions that follow from a Landau Ginzburg analy- 
sis to be discussed elsewhere. We obtain the curve 
X _1 (T) (which is plotted in the inset of Figj|), that can 
be approximated, close to the transition, with the ex- 
pression x' 1 — ~^bi wnere t i s the reduced temperature 
t = T — T c and a and b are real positive parameters. 
More precisely (d)(U,T) contains a singular part (aris- 
ing from the liquid gas analogy close to the transition) 
and a regular part. The regular part gives the constant 
background (b/a) and the singular part produces the di- 
vergent susceptibility (a/t) which is linear in the inverse 
reduced temperature. T c and U c are fitting parameters, 
which give the location of the second order critical point 
where the first order lines terminate. A least squares fit 
to x gives our estimate for T c = 0.026 ± 0.003 and us- 
ing the liquid gas analogy we find that U c = 2.38 ± 0.02. 
Thus, we have obtained estimates for the position of the 
finite temperature second order critical point in the sin- 

tband Hubbard model using the QMC method ]2C| ] 
. It is important to realize now that the lowest tem- 
perature set of data for the double occupation in Fig jl] at 
T =1/40 = 0.025 appear to be slightly below T c . These 
results thus indicate a discontinuous jump in the double 
occupation. 

A more general consequence of the presence of the fi- 
nite temperature second order critical point, is the ex- 
istence of a region at T lower than T c , where two dif- 
ferent solutions of the mean field equations exist. This 
coexistence of solutions, as usual, results in first order 
transition lines. We have, therefore, searched for differ- 
ent converged solutions well below T c = 0.026 and per- 
formed calculations at T = 1/51.2 and T = 1/64 for 
several values of U close to U c . Indeed, we obtained two 
different and fully converged solutions at U = 2.4 in the 
first case, and at U = 2.4 and U = 2.5 in the latter. For 
the U = 2.4 and T = 1/64 solutions, we have taken spe- 
cial care to rule out the possibility of systematic errors, 
by performing a large number of further iterations after 
self-consistency was achieved and also by increasing the 
number of sweeps to up to 300, 000 to minimize statistical 
fluctuations. In Figj^, we display two coexistent solutions 



at U = 2.4 and T = 1/64, along with the results for the 
imaginary part of the Green's functions on the real fre- 
quency axis using the maximal entropy method p2] for 
the analytic continuation. These results, obtained within 
a narrow range of interaction U close to U c and at tem- 
peratures below our estimate for T c , further confirm the 
consistency of our numerical results. 

In conclusion, we presented the results of a careful nu- 
merical QMC study of the finite temperature Mott tran- 
sition in the paramagnetic phase of the Hubbard model 
in infinite dimensions. We identified the singular behav- 
ior of a susceptibility associated with the finite temper- 
ature critical point and good estimates for U c and T c 
were obtained. This should be contrasted with a previous 
QMC study of the same problem Jll[] which reached the 
conclusion that there is no signature of a finite temper- 
ature second order phase transition in the temperature 
region accessible to QMC studies. Furthermore they pre- 
sented Tbound — 0.01429-D as a rigorous numerical upper 
bound on T c . Here we have shown using two different 
approaches, i.e., a high temperature susceptibility cal- 
culation and an explicit demonstration of coexistence at 
low temperatures, that a finite temperature critical point 
exists and that the bound obtained in |flj is incorrect. 

In earlier publications we asserted that qualitative fea- 
tures of the paramagnetic metal to paramagnetic insula- 
tor transition observed in V2O3, are well described by the 
DMFT of the single band Hubbard model on a frustrated 
lattice treated using the iterated perturbation theory. In 
this work we show that those conclusions are also valid 
when a more accurate technique, such as QMC, is used 
to solve the dynamical mean field equations. In partic- 
ular, the anomalous temperature dependence of physical 
properties (double occupancy, resistivity etc.) and the 
crossovers discussed in Ref. ||, which were tied to the 
proximity to a second order finite temperature metal to 
insulator transition in the model. 

While the qualitative agreement between the DMFT 
and the experiments is remarkable, there are quantita- 
tive disagreements. The measured plasma frequency is 
a factor of two larger than the calculated one (using ex- 
act diagonalization and DMFT), if we take the LDA esti- 
mate of the bandwidth D w 0.5eV p4f . This was already 
noticed in Refs. H^3|. Using the same LDA based esti- 
mates for the model parameters, the critical temperature 
that we compute (with DMFT and QMC) is also approx- 
imately a factor of two smaller than the observed T c in 

v 2 o 3 §. 

These results suggest the necessity for the inclusion of 
additional features like orbital degeneracy, ligand bands 
p5| and the electron phonon interaction 26 in the model 
to be able to make accurate predictions for physical quan- 
tities in this system. The extension of the QMC algo- 
rithm to orbitally degenerate models large d models pro- 
posed by one of us |27| should be useful for this purpose. 
Indeed, less extensive studies than those carried out al- 
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ready show that orbital degeneracy affects the critical 
temperature [ p7[ . 
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FIG. 1. The double occupation (d) as a function of U for 
several values of T= 1/20, 1/25, 1/28 1/32, 1/35, 1/40 (from 
bottom to top). The inset shows the inverse susceptibility 
X -1 as a function of T. The line is a least squares fit using 
the expression in the text. The intercept with the T— axis 
gives our estimate for T c = 0.026 
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FIG. 2. Two converged solutions, one metallic-like (solid 
circles) and one insulating- like (open circles), of the DMFT 
equations for the same value of the interaction U = 2.4 and 
T — 1/64. In the inset we show their corresponding density 
of states. 
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